// project home: cp3-cg.googlecode.com

#include "myvectormath.h"
#include <iostream>
using namespace std;

vector<Double_t> matrixvectorproduct(vector<Double_t> *A, vector<Int_t> *JA, vector<Int_t> *IA, vector<Double_t> *x, bool debug){
  
//   bool debug=true;
  
  Int_t dim = x->size();
  if (debug){cout << "dimension: " << dim << endl;}
  vector<Double_t> b(dim,0); //Ausgabevektor
  
  for (int j=0;j<dim;j++){
    if (debug){cout << "j=" << j << ";" << endl;}
    for (int i=IA->at(j);i<IA->at(j+1);i++){
      if (debug){cout << "i=" << i << ";" << endl;}
      b.at(j)=b.at(j)+A->at(i)*x->at(JA->at(i));
    }
  }
  
  return b;
}



Double_t skalarp(vector<Double_t> *x, vector<Double_t> *y){
  Double_t product=0;
  for (int i=0; (unsigned)i<x->size();i++){
    product=product+x->at(i) * y->at(i);
  }
  return product;
}


vector<Double_t> vectorsub(vector<Double_t> *x, vector<Double_t> *y){
  vector<Double_t> diff(x->size(),0);
  for (int i=0; (unsigned)i<x->size();i++){
    diff.at(i)=x->at(i) - y->at(i);
  }
  return diff;
}


vector<Double_t> vectoradd(vector<Double_t> *x, vector<Double_t> *y){
  vector<Double_t> sum(x->size(),0);
  for (int i=0; (unsigned)i<x->size();i++){
    sum.at(i)=x->at(i) + y->at(i);
  }
  return sum;
}


vector<Double_t> VectorMultiS(vector<Double_t> *x, Double_t alpha){
  vector<Double_t> prod(x->size(),0);
  for (int i=0; (unsigned)i<x->size();i++){
    prod.at(i)=x->at(i) * alpha;
  }
  return prod;
}





